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We use molecular dynamics simulations to study the exchange kinetics of water molecules at a model metal 
electrode surface - exchange between water molecules in the bulk liquid and water molecules bound to the 
metal. This process is a rare event, with a mean residence time of a bound water of about 40 ns for the 
model we consider. With analysis borrowed from the techniques of rare-event sampling, we show how this 
exchange or desorption is controlled by (1) reorganization of the hydrogen bond network within the adlayer of 
bound water molecules, and by (2) interfacial density fluctuations of the bulk liquid adjacent to the adlayer. 
We define collective coordinates that describe the desorption mechanism. Spatial and temporal correlations 
associated with a single event extend over nanometers and tens of picoseconds. 


I. INTRODUCTION 

Hydrated metal interfaces are standard platforms 
for heterogenious catalysis.® Understanding their micro¬ 
scopic kinetic processes is important for the design of 
efficient and clean catalytic systems.® Currently, this un¬ 
derstanding is limited by too little knowledge of the 
molecular-level structure and dynamics at interfaces. 
Most experimental work has come from macroscopic 
electrochemical measurements where inferring molecular- 
level details is challengingor surface science measure¬ 
ments at ultrahigh vacuum conditions, where processes 
relevant to condensed phase behavior are absent.®^ Most 
theoretical work has utilized quantum chemical meth¬ 
ods designed to produce highly accurate representations 
of one (or a few) configurations, like that reviewed in 
Refs. [2]and|7l That type of study neglects fluctuations, 
relaxation and solvation. Here, we adopt a complimen¬ 
tary perspective, following in the footsteps of Refs. 
employing classical molecular dynamics simulations to 
study large-length scale correlations and long-time scale 
dynamics. 

We use these methods to study the dynamics of ex¬ 
change between water molecules in the bulk liquid with 
water molecules tightly bound to a platinum electrode. 
We show how the bound water molecules organize into 
an adlayer with long-lived disordered spatial patterns ex¬ 
tending over large distances. Similarly, we show how the 
bulk water molecules adjacent to the monolayer form a 
soft liquid interface that can accommodate long wave¬ 
length distortions in the local interfacial density. The 
exchange between the bulk and the bound waters involves 
the dynamics of both the adlayer and the soft liquid inter¬ 
face. The exchange process is rare, occurring with a char¬ 
acteristic timescale on the order of tens of nanoseconds. 
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It is collective, proceeding through a reaction coordinate 
that involves reorganization of many adlayer molecules 
in concert with fluctuations of the soft liquid interface. 

This work highlights that in the vicinity of an aqueous 
electrode interface, even the simplest kinetic processes 
become highly collective. It is in line with our previous 
work that showed how an aqueous ele ctrode interface ex¬ 
hibits an array of emergent behaviorsIn particular, 
large binding energies pin water molecules to the surface 
resulting in a water adlayer in which hydrogen bonding 
is passivated almost entirely within the plane of the elec¬ 
trode. This in-plane hydrogen bonding greatly reduces 
interactions between the adlayer and surrounding liquid, 
producing a composite water-metal surface that is hy¬ 
drophobic. The composite surface attra cts defects from 
the bulk li^id like excess protons ,1^^^ and weakly sol¬ 
vated ions,^ and it facilitates large density fluctuations 
in the adjacent bulk liquid.^ Further, when surface peri¬ 
odicities and local coordinations of the electrode surface 
are incommensurate with favorable hydrogen bonding ge¬ 
ometries, hydrogen-bonding patterns of the adlayer are 
frustrated, leading to structures with localized defects 
and heterogeneous dynamics 

Consequently, the dynamics at this interface is gov¬ 
erned by a hierarchy of timescales, spanning picoseconds 
to tens of nanoseconds. The shortest of these timescales 
characterize molecular fluctuations in the bulk liquid, 
the longest characterize relaxation of the adsorbed water 
monolayer These fast and slow degrees of freedom are 
coupled through fluctuations of the bulk-monolayer in¬ 
terface that itself relaxes over intermediate timescales 
This complexity affects the exchange of a water molecule 
between the adsorbed monolayer and the bulk, and as 
such, it also affects the accessibility of reactive sites on 
a hydrated electrode. Recognizing this fact, Hansen et 
al introduce phenomenological energy barriers in their 
treatment of oxygen evolution.^^^ These barriers are at¬ 
tributed to entropic effects related to surface water reor¬ 
ganization, highlighting the role of slow interfacial water 
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dynamics. Related effects may be important in gas phase 
catalysis as well because a similar monolayer of water 
forms on platinum surfaces exposed to ambient air. Slow 
reorganization of this monolayer has been conjectured 
to account for the low sticking coeffi cient of gaseous ad¬ 
sorbates such as carbon dioxidedecreasing catalytic 
efficiency. 

Figure illustrates the basic steps in the first half of 
a typical exchange process. This first half is the desorp¬ 
tion of a water molecule from the electrode surface as 
found in the simulations we detail in the next sections of 
this paper. The three panels of the figure are each sepa¬ 
rated by approximately 1 ps, with a tagged molecule ini¬ 
tially adsorbed to the metal, proceeding through its tran¬ 
sition state where it leaves behind a surface vacancy, be¬ 
fore committing to bulk liquid where the vacancy either 
diffuses away or becomes annihilated by an adsorption 
event. The orientational relaxation of water molecules 
within the adlayer occurs roughly 10^ times more slowly 
then for those within the bulk liquid, and nearly 10^ times 
more slowly than at the bulk liquid interface. Not sur¬ 
prisingly, a detailed characterization of molecular desorp¬ 
tion requires the quantification of local surface reorgani¬ 
zation ability. 

Figure also renders the liquid’s instantaneous inter¬ 
face adjacent to the composite electrode-water surface. 
This soft-liquid interface is identified as in Ref. [22l We 
see below that it is important in distinguishing the ad¬ 
sorbed and desorbed states with reference to the posi¬ 
tion of the instantaneous liquid phase boundary. Only 
with this frame of reference does it become clear that 
desorption is an instantonic event with few recrossings. 
Otherwise the nature of the process is obscured by slow 
capillary-wave-like dynamics of the interface. With the 
adsorbed and desorbed states or basins for a tagged wa¬ 
ter molecule thus defined, we use tools of transition path 
samplin^^ to consider the mechanism of desorption. In 
particular, we determine a collective reaction coordinate 
that describes the ensemble of transition states. Finally, 
we consider the dynamics involved in the time-reverse 
process, i.e., the adsorption of a molecule from the bulk 
onto the electrode. The Appendix provides additional 
results and details on our methodology. 


II. BASIN DEFINITIONS AND TIMESCALES 

In order for an adsorbed molecule to enter the bulk liq¬ 
uid it must incorporate itself into the established hydro¬ 
gen bonding environment of the liquid. Such a process 
is initiated through the scavenging of broken hydrogen 
bonds that arise spontaneously but fleetingly at a liquid 
water interface. This process is thus modulated by spa¬ 
tial fluctuations in the position of the bulk liquid inter¬ 
face that, due to the hydrophobic nature of the extended 
water-metal interface,^ are large and driven by collective 
displacements of weakly-bound molecules. Consequently 
if the position of the desorbing molecule is described rel- 



FIG. 1. (left) Characteristic snapshot the adsorbed state. The 
tagged molecule, highlighted in yellow, sits atop a platinum 
atom (grey), below the instantaneous interface (blue), and hy¬ 
drogen bonds to other waters (red and white) predominately 
within the adlayer plane, (center) Characteristic snapshot 
of a member of the transition state ensemble. A rare density 
fluctuation in the liquid above, rendered with the deformed in¬ 
stantaneous interface, allows for the tagged molecule to come 
off of the electrode, leaving a vacancy on the surface, shown in 
yellow, (right) Characteristic snapshot of the desorbed state. 
The tagged molecule is above the instantaneous interface and 
the vacancy diffuses away. Snapshots are separated by 1 ps. 
A movie of the same trajectory can be streamed online at 
https://youtu.be/MVjUBTBrmfU. 


ative to a static frame of reference, the fluctuations in 
the soft liquid interface are likely to obscure meaningful 
details related to the formation of the initial hydrogen 
bonds between a desorbing molecule and the bulk. In or¬ 
der to account for these fluctuations explicitly, we adopt 
a reduced description of the instantaneous density field 
that incorporates these fluctuations. 

The instantaneous liquid interface, s, is identified im- 
plicitlyP^ 

p{s) = Ptl2 , (1) 

where pi is the bulk liquid density, and p(r) is the molec¬ 
ular density coarse grained over a length scale i.e., 

N 

= ( 2 ) 

i 

Here, is the position of the oxygen of the Rh wa¬ 
ter molecule, there are N molecules in total, and 0(r) 
is chosen to be a spherically symmetric and normalized 
Gaussian function, truncated and shifted to become zero 
continuously beyond at |r| = 3^. We take ^ = 2.5A. 

The position of the Rh molecule relative to the inter¬ 
face is 

ai = n(Sri) • i^i - SrJ , (3) 

where n(s) is the unit vector normal to the interface at 
s pointing away from the bulk, and Sr is the point on s 
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FIG. 2. Distributions and timescales for desorption, (a) Po¬ 
tential of mean force for moving a water molecule relative to 
the instantaneous interface. The red and black dashed lines 
denote the maximum and minimum values for defining basins 
of adsorbed and desorbed states. In P{a) is shifted by a con¬ 
stant so as to make its minimum zero, (inset) Characteristic 
snapshot of an adsorbed molecule and the instantaneous liq¬ 
uid interface. Flux correlation function, blue solid line, as 

defined in Eq. The dashed red line is a linear fit to C{t) 
whose proportionality constant yields the rate for desorption, 
(c) Distribution of instanton times for the desorption event, 
blue, as defined in Eq. and a reference exponential decay 
with characteristic time At = 3ps. 


closest to r. See inset to Fig. is a scalar that can 

be positive or negative. We employ the convention that 
omits the tagged molecule from the density used to de¬ 
fine the instantaneous interface. The alternative, includ¬ 
ing the tagged molecule in the density, gives behaviors 
qualitatively similar to those shown below. 

With an umbrella sampling procedure developed previ¬ 
ously ,1^ we calculate P{a) = ((5(a — a^)), where the angle 
brackets, (•••), denote equilibrium average. The func¬ 
tion — In P(a) is the associated potential of mean force in 
units of Boltzmann’s constant times temperature, /cbT. 
This function is shown in Fig.j^a) for water adjacent to 
the 100 Pt surface. It has three notable features: (1) For 
a < 0, there is a narrow minima reflective of the bound 
adlayer. (2) Between the bound adlayer and the bulk liq¬ 


uid, there is a free energy barrier peaked at a ^ 0 with a 
height relative to the surface layer of 4 /cbT. (3) The bulk 
liquid, a > 2 A, has a more diffuse basin than that of the 
bound layer. Similar qualitative features appear in the 
analogous free energy as a function of absolute position 
shown in the Appendix. 

In view of these features, we define the following dy¬ 
namic indicator functions for a tagged molecule, i, be¬ 
longing to the adsorbed state, hi, or desorbed state, gi. 
These are defined as 

hi = 1, if ai < -0.7 A (4) 

= 0, else 

and 

5i = l, if ai>0.3A (5) 

= 0, else 

and are indicated in Fig.j^a). The time-correlation func¬ 
tion between the two 

(«) 

is shown in Fig. [^b). We use standard notation for the 
time correlation function, i.e., the equilibrium average is 
over initial conditions, and hi{t) is the desorbed state in¬ 
dicator function evaluated for the configuration at time 
t. The correlation function is independent of the particle 
label, i, because all the water molecules are indistinguish¬ 
able at equilibrium. Later, we generalize to nonequilib¬ 
rium cases where the particle-label index is important. 

This correlation function is estimated with straight¬ 
forward molecular dynamics simulations by harvesting 
2000 trajectories bridging adsorbed and desorbed basins 
(see Appendix for details). The short-time behavior of 
C{t) is indicative of transient barrier crossings and re¬ 
crossings and the longer-time behavior exhibits the linear 
growth expected for a rate process .1^ The time deriva¬ 
tive, dC{t)/dt, reaches its plateau value after a transient 
time, and this transient time is approximately the mean 
instanton time, (At), where a formula for the instanton 
time of a specific trajectory. At, is given below. With the 
basins as defined above (Eqs.0and[^, this time is about 
3 ps. Far beyond that time^.e. t ^ (At)), the slope 
determines the rate constant, i.e., C{t) ^ k^t. The best 
fit to the data in Fig. [^b) yields /cd ~ 2.6 x 10“^ ns“^. 
In other words, given the surface density of water on the 
electrode, one desorption event occurs roughly every 3 ns 
per nm^. 

The rate for desorption is independent of the basin def¬ 
initions, but the typical time to reach the plateau, (At), 
is sensitive to these definitions. This sensitivity arises 
from the possibility of recrossing events as the molecule 
lingers near the top of a free energy barrier. As a quanti¬ 
tative measure of the time it takes to cross the barrier, for 
each single trajectory of duration tobs which the tagged 
molecule i crosses from adsorbed state to desorbed state. 
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we compute the instanton time, 

rtohs 

At= dt [l-gi{t)] . (7) 

Jo 

From an ensemble of those trajectories, we obtain the dis¬ 
tribution of instanton times, which is plotted in Fig.j^c). 
The distribution is roughly exponential, consistent with 
different desorption events being uncorrelated. The time 
constant for the exponential is about 3 ps. 

This average value, (At), is in line with the typical time 
for a water molecule to diffuse lA in the homogeneous 
liquid, which is about 5 ps. On the other hand, a much 
larger instanton time is found when using definitions of 
hi and gi based upon the separation from the metal sur¬ 
face or mean interface, rather than the separation from 
the instantaneous interface. That alternative separation 
is given by the component of the vector perpen¬ 
dicular to the plane of the static surfaces. In that case, 
replacing ai with Zi in Eqs.[^and[^ we find that the av¬ 
erage instanton time increases to 20 ps. (See Appendix 
B.) This larger time indicates the increased likelihood 
for recrossing events along this coordinate. Indeed this 
longer timescale is commensurate with the average time 
for height fluctuations of the instantaneous inter face, 
suggesting that the physical origin for the additional re¬ 
crossing events lie in the relaxation of these interfacial 
fluctuations. 


III. ANALYZING THE TRANSITION STATE 
ENSEMBLE 

The variable - the separation of a water molecule 
from the instantaneous interface - serves as a distinguish¬ 
ing order parameter for adsorbed and desorbed states, 
but it is not sufficient to characterize pathways between 
adsorbed and desorbed states. This fact is demonstrated 
now through analysis of the transition state ensemble 
(TSE). 

A. Committor and definition of TSE 

Members of the TSE can be harvested from the tra¬ 
jectories used to compute the correlation function, C{t). 
These are the configurations from which short trajec¬ 
tories commit to either reactant or product with equal 
likelihood.l^ In particular, for trajectories that begin at 
coordinate x = {ri,...,rjv} with randomly chosen ve¬ 
locities and propagate for a time tobs, the commitment 
probability, p(x, tobs): is the fraction of trajectories that 
end up in the desorbed basin. Here, the trajectory time 
is larger than the transient time but not so large to per¬ 
mit multiple transitions, i.e. (At) < tobs ^ l/^d, and 
accordingly, the TSE is the subset of configurations with 
p(x, tobs) = 1/2. To estimate p(x, tobs), we use tobs = 10 
ps and average over 10 separate realizations of the ran¬ 
dom velocities taken from the Maxwell-Boltzmann dis¬ 


tribution. If this estimate falls within a 90% confidence 
interval of 0.5, we say x is a member of the TSE. 

The distribution of these configurations, Ptse(x), is 
difficult to interpret because x is a high-dimensional 
variable. To make headway, we project onto a low¬ 
dimensional representation. The notation we use for the 
resulting marginal distributions is 

-Ptse(- 4) = y dxPTSE(x) <5 (i(x) - , (8) 

where A(x) is the dynamical variable, and A is its ob¬ 
served value. 

Eigure [^a) refers to the TSE projected onto the co¬ 
ordinate ai, i.e., the separation of the desorbed molecule 
from the instantaneous bulk liquid interface (see Eq. [^. 



FIG. 3. (a) The probability distribution for the relative dis¬ 

tance to the instantaneous interface, for a tagged molecule 
within the transition state ensemble. In grey is the equilib¬ 
rium distribution, (b) The probability distributions for the 
absolute position for a tagged molecule within the transition 
state ensemble. In grey is the equilibrium distribution, (c) 
Three free energy profiles relative to the instantaneous in¬ 
terfaces (in units of ksT): for moving a water molecule re¬ 
versibly (black); for moving a water molecules either from the 
surface or from the bulk, reversibly with respect to the bulk 
waters, but irreversibly with respect to the other surface ad- 
layer molecules, which are fixed in one typical configuration 
(blue or red, respectively). 
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A/ZcbT X = F/kBT, G/ksT 


FIG. 4. (a) Probability distribution of surface reorganization energies, A for a given surface configuration q. The dashed line 
is a Gaussian fit. (b) Probability distributions of the work associated with pulling an adsorbed molecule off of the electrode, 
with surface configuration q, and under different constraints. Dashed lines are Gaussian fits, (c) Gharacteristic snapshot of 
adsorbed waters, with metal atoms colored by their value of surface reorganization energy, A. 


Figure [^b) refers to the TSE projected onto the co¬ 
ordinate Zi, i.e., the separation from the average posi¬ 
tion of the adlayer. Both distributions are peaked in 
the interfacial region between the bound adlayer and the 
bulk liquid. The relative widths of the distribution are 
{\6a\)/{a) ^ 0.2 and {\5z\)/{z) « 0.35. By this met¬ 
ric, reference to the instantaneous interface (i.e., the dis¬ 
placement a^), is a better descriptor of the TSE than 
reference to the mean interface (i.e., the displacement 
Zi). However, in both cases, Eigs. |^a) and (b) show, 
the transition states still occur over a significant range 
of separations. These distribution widths indicate that 
the desorption reaction coordinate has important com¬ 
ponents orthogonal to both of these distances. 

B. Role of adlayer reorganization 

Important variables missing from (or zi) involve con¬ 
figurations of the water molecules bound to the metal 
surface. We use the symbol q to stand for the collection 
of coordinates that together with specify the config¬ 
uration of that adlayer. q is significant because each 
desorption event is accompanied, at least transiently, by 
the creation and relaxation of a surface vacancy. Effects 
of this reorganization can be seen by juxtaposing the re¬ 
versible work function, — In P(a), with two nonequilib¬ 
rium counterparts that depend upon q. 

The first counterpart is 

G'i,q(a,ao) = da'(/i(a'))q. (9) 

Here, fi{a) is the instantaneous force on evaluated at 
ai = a [i.e., fi{a) = -{dU/dai)\a,=a, where U = P(x) 
is the total potential energy of the system], and (• • • )q 
stands for equilibrium average with q fixed. Gi,q(a,ao) 
is therefore the reversible work (in units of /cpT) to move 
ai from a reference position, ao, to a with the rest of 
the adlayer fixed in configuration q. Averaging this con¬ 
strained work over the equilibrium distribution of q yields 


— In P(a) to within an additive constant. Without that 
averaging, the constrained free energy function is gen¬ 
erally different from the reversible work function. Two 
typical examples of the constrained free energy function 
are illustrated in Eig. [^c) with ao = — 1 A (the dashed 
blue line) and with uq = 1 A (the dashed red line). The 
respective q’s are taken from the equilibrium ensembles 
of configurations with set at a value where = a^. 
Comparison with — InP(a) (the black line) shows that 
adlayer reorganization is indeed a significant and irre¬ 
versible effect. 

Some of that reorganization occurs on relatively short 
time scales, which brings us to the second nonequilibrium 
counterpart, 

Fi,q{a,ao;t) = j da'(/i(a'))q,t • (10) 

Here, (• • • )q,t stands for the time average over a trajec¬ 
tory of length t, with q as the initial configuration for 
the rest of the adlayer. q, in this case is not constrained 
in the subsequent dynamics. Pi,q(a, ao; t) is then the av¬ 
erage work done to move a^ from ao to a in a time t 
starting with the rest of the adlayer in configuration q. 
As t ^ oo, this time-dependent work function becomes 

— InP(a) to within an additive constant. 

On the way to that limiting behavior, there are two re¬ 
laxation regimes. The first coincides with reorganization 
directly spurred by the desorption event. It takes place 
on times scales between 10 and 100 ps. The second coin¬ 
cides with equilibration of the adlayer’s transient hetero¬ 
geneity. It takes place on time scales of 1 to 100 ns. Thus, 
P^^q(a, ao; t) for 1 ns < t < 10 ns gives a measure of relax¬ 
ation for one realization of the transient heterogeneity, 
and as such, it has a distribution of values characteristic 
of the distribution of transient heterogeneity. 

In general, for fixed a and Uq, G^^q and P^^q are differ¬ 
ent due to the different constraints on the relaxation of 
surface bound molecules. When G^^q and P^^q are com¬ 
puted with ao near the equilibrium position of adsorbed 
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state and a in the bulk, their deviation signals the addi¬ 
tional local free energy that must be dissipated in order 
to accommodate a molecule leaving the surface. In this 
way, by locating tagged regions where both measures are 
equivalent, we can isolate regions on the surface where re¬ 
organization is facile and admitting to desorption events. 
Specifically, we can define the deviation between these 
two work functions for moving a tagged molecule between 
the surface and bulk for a given initial surface configura¬ 
tion 



X = G-F. (11) 

We adopt the shorthand, G = Gi,q(a,ao) and F = 
F^^q(a,ao) where each is evaluated with ao = —1 and 
a = 3. To determine the reorganization energy dissi¬ 
pated through local surface relaxation, our results are 
not sensitive to these precise values as long as uq is near 
the maximum of P{a) and a is within the large a plateau. 
A histogram of A, denoted Pq(A) to emphasize that each 
value of A depends on the initial surface configuration, 
is shown in Fig. for the 100 surface. This is deter¬ 
mined from the histograms of G and F shown in Fig.|^. 
Within our limited statistics each distribution is Gaus¬ 
sian. The peak value of G is larger than F by 7 /cbT and 
its distribution is wider by a factor of 1.5. The larger av¬ 
erage value of G over F is a consequence of the additional 
constraint on the surface. The variance in A is not the 
sum of the variances of G and F, implying some degree 
of correlation, however this correlation is small, with a 
correlation coefficient of 0.2. Distributions of A, G and F 
for the 111 surface are shown in the Appendix. 

For the 100 surface, we find that on average the local 
number of persistent hydrogen bonds made adjacent to 
a tagged site increases with increasing A. We define a 
persistent hydrogen bond as an oxygen oxygen distance 
of less than | Too |= 3.5A and a OH bond angle rela¬ 
tive to Too less than 30^ with positions coarse-grained 
over 1 ps. This is a reflection of the stable monolayer 
surface structure, which is commensurate with the 4- 
fold coordination of the exposed lattice and water’s pre¬ 
ferred hydrogen bonding patterns. Defects that disrupt 
this four-coordinated network allow for increased mobil¬ 
ity of the water molecules orientation,E^and consequently 
these sites have a lower reorganization energy. However, 
there is not a perfect correlation between persistent hy¬ 
drogen bonds and reorganization energy, the correlation 
coefficient is 0.5. This low correlation is because the reor¬ 
ganization energy is not a local quantity, and rather can 
reflect more subtle strains in the hydrogen bonding net¬ 
work that add up over neighboring molecules and allow 
for facile relaxation. Hydrogen bonds drawn in Fig. [^c) 
are the persistent hydrogen bonds discussed above. The 
trend relating hydrogen bonding and A can be observed 
from this snapshot, though also shown are exceptions 
such as molecules with only 3 persistent hydrogen bonds 
that can still have relatively high values of A. While 
many-bodied in nature, as illustrated in Fig. [^c) the 
spatial distribution of reorganization energies are never¬ 


FIG. 5. Commitor distribution function averaged over an 
ensemble constrained to A ~ 0 for the (a) 100 and (b) 111 
surfaces. The red dashed lines are Bernoulli distributions for 
11 trials. 


theless made up of uncorrelated domains that vary over 
the surface with negligible correlation lengths, similar to 
that found for the induced charge distributions on the 
metal centers .1^ 

As an ultimate test of the relevance of this surface co¬ 
ordinate for the reaction dynamics, we have calculated 
the distribution of commitment probabilities over an en¬ 
semble initialized with A ~ 0. This type of test, first 
proposed by Geissler, et. aP^ distinguishes potential re¬ 
action coordinates that well characterize the TSE only on 
average due to a near equivalency of populations of con¬ 
figurations committed solely to the A or B basins. Figure 

plots the distribution of commitor values for fixed A, 
Px{Pb) with a reference to a binomial distribution with 11 
trials, corresponding to the mean number of trials used in 
the evaluation of the commitment probability. The peak 
near 0.5 confirms that these constraints impose a prox¬ 
imity to the true seperatrix, and their width is consistent 
with a finite number of trials. Figure plots these dis¬ 
tributions for both the 100 and 111 surfaces of platinum, 
illustrating the degree of generality of this reaction co¬ 
ordinate. Details on the calculations for the 111 surface 
are provided in the Appendix. 


IV. CORRELATIONS BETWEEN EXCHANGE PAIR 

Most of this work has focused on the microscopic 
events leading a water molecule to desorb from the elec¬ 
trode surface. By microscopic reversibility, the pathways 
for the adsorption of a water molecule from the bulk onto 
an electrode must be the time-reversed as those for the 
desorption of a water molecule from the electrode into the 
bulk. Moreover, from the law of mass action, the rate of 
adsorption must be equal to the rate for desorption times 
the free energy difference between the two states. There¬ 
fore, it would seem that by studying the forward process 
of desorption, all relevant information of exchange can be 
extracted. Though this is mostly true, neither of the two 
previous statements prohibits the existence of dynamic 
correlations between the exchange pair. 
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FIG. 6. Characteristic time series for vacancy diffusion during an exchange process. The exchange pair are highlighted 
in yellow as is the surface vacancy. Snapshots are approximately 5 ps apart. A movie of this trajectory found here: 
https://youtu.be/P8LfvycgdpI. 


The extent of correlations between the exchange pair 
is evident in Fig. which depicts a characteristic time 
series of the diffusion of the vacancy left behind after a 
desorption event. The kth vacancy is created each time 
a desorption event occurs, at the position r, the location 
of the closest platinum atom to the desorbing molecule, 
is denoted n/e(r,t) and carries a numerical value of 1. 
Similarly, a vacancy is annihilated each time an adsorp¬ 
tion event occurs, which changes the numerical value of 
n/c(r,t) to 0. A typical vacancy is highlighted in yellow 
in Fig. and is shown to diffuse through the surface over 
the course of 30 ps before being annihilated by an adsorb¬ 
ing molecule. The diffusion of the vacancy requires the 
hydrogen bonding network of the adlayer to reorganize, 
which is does following a process similar to that described 
by Ref. |28]in low water coverage ab initio simulations on 
the palladium 111 surface, whereby a hydrogen bonded 
pair is anchored by one water on a metal center and ro¬ 
tates its hydrogen bonded partner that is weakly bound 
to the surface between two adjacent sites. 

The lifetime of a surface vacancy determines the con¬ 
ditional waiting time for an adsorption event, given a 
desorption event just occurred, in a way that is local and 
not obscured in the limit of large surfaces. We define the 
tex, as the lifetime for the kth vacancy. 


t 


ex 


dtn/e(r,t). 


( 12 ) 


Figurej^a) shows the distribution of lifetimes in logscale, 
which is strongly peaked at, (tex) =60 ps. This time 
is much less then the mean time for adsorption, which 
means most adsorption events result from correlated sur¬ 
face configurations that transport localized defects from 
one place on the surface to another. 

The pair correlations for the location of the creation 
and annihilation of the surface vacancy mirror the va¬ 
cancy lifetimes. Specifically, we define a pair correlation 
function for the location of the vacancy creation and it 
subsequent annihilation. 


/)„„(|r-r'|) 


{nk{r,t)nk{r',t + te^)) 
(nfc(r',i + iex)) 


(13) 


which due to rotational invariance depends only on the 
magnitude of r — r' , and due to time translational invari¬ 
ance, depends only on tex- This correlation function are 
shown in Fig.[^b). While it is most probable that the va¬ 
cancy does not diffuse, the decay in Fig.j^b) is slow, and 
there is significant probability that vacancies diffuse over 
nanometer lengths. This function can be well fit with a 
biexponential decay with characteristic lengths of 0.5A 
and 4.5A. The long waiting time and large distance cor¬ 
relations can be understood from what has already been 
found regarding the kinetics of desorption. Namely, ad¬ 
sorption occurs due to a localized defect in the hydrogen 
bond network of the surface waters coming into contact 
with a rare, fleeting fluctuation of the interfacial density. 



FIG. 7. (a) Probability distribution of lifetimes for a surface 
vacancy. The mean lifetime is (tex) = 60 ps. (b) Probability 
distribution for the distance between the surface sites for va¬ 
cancy creation and the following site of vacancy annihilation. 
The dashed red line represents an exponential distribution of 
exchange times, which on the the logarithmic scale has a char¬ 
acteristic single peak at its mean value. Gharacteristic decay 
lengths are 0.5A and 4.5A. 
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Appendix A: Methods 

The model we adopt to describe the interface utilizes a 
standard classical force field to compute the intermolec- 
ular interactions between water molecules, namely the 
SPC/E potential,!^ and a semi-empirical description of 
an ideally polarizable metal surface held at constant elec- 
trical potential developed by Seipman and Sprlk.^^SEH 
The atomic geometry of the electrode, along with the 
details of the water-metal interactions are consistent 
with that of crystalline platinum. The water-metal in¬ 
teraction is based on simple quantum-chemical calcu¬ 
lations and reproduces the large binding energy and 
ground state geometry of a single water molecule.l^ Clas¬ 
sical contri butio ns to the potential of zero charg^^^ and 
capacitanc d^^ l ^^ l have been shown to be in good agree¬ 
ment with more detailed models. 

The system dimensions and particle numbers are the 
same as in previous work. Specifically, a 3.5nm by 3.5nm 
by 5nm simulation box is filled with nearly 2000 water 
molecules and close to 1000 platinum atoms in 2 slabs, 
each three rows wide, resulting in a mean bulk water den¬ 
sity in the center of the slab of 1 g/cc. Equilibration is 
difficult due to the slow surface water dynamics, so an 
ensemble of initial conditions is generated by annealing 
the system at 400 K over 1 ns, and cooled to 298 K over 
another ns. Production runs used to generate the en¬ 
semble of reactive trajectories are 20 ns, which for each 
20 surfaces yield approximately 2000 desorption events. 
The integrator uses a langevin thermostat, with charac¬ 
teristic time constant of Ips, and molecular constraints 
imposed with SETTLE.Long ranged electrostatics are 
accounted for using the appropriate Ewald summation 
for mixed point and Gaussian charge distributions,l^and 
the two-dimensional periodicity is approximated with the 
standard slab correction,!^ which elongates the simula¬ 
tion cell by a factor of 5. We use the same methodology 
for both 100 and 111 surfaces. Galculations are done with 
a modified version of the LAMMPS code.!^ 

Umbrella sampling with MBARP^ is used to estimate 


FIG. 8. Potential of mean force for moving a water molecule 
relative to the absolute distance from the electrode surface, 
(inset) Flux correlation functions for z basin definitions for 
the 100, red, and 111, black surfaces. In blue is a line with 
unit slope. 


free energies in Figs. [2][4| For free energy calculations 
with quenched surface configurations, production runs of 
1 ns are sufficient to converge local distributions, while 
for calculations with evolving surface configurations pro¬ 
duction runs of 10 ns are required. In each case, 12 win¬ 
dows are equally spaced throughout the a coordinate and 
harmonic biasing potentials with spring constants of 8.0 
kcal/mol A are used. 

Appendix B: Cartesian coordinate basin definition 

An alternative means for defining the basin definitions 
in Eq. and would be to use the absolute position 
away from the electrode, rather then the position rela¬ 
tive to the instantaneous liquid interface. To construct 
the analogous indicator functions it is first useful to cal¬ 
culate the free energy for moving a molecule in along this 
coordinate. This free energy is easily calculated by tak¬ 
ing the negative logarithm of the average water density 
as a function of 2 ;, or p{z). 

Eigure plots the free energy as a function of z for 
both the 100 and 111 surface. Similar to the results 
in Eig. [^a), both curves exhibit a narrow basin near 
the electrode, z < lA, and a large barrier between this 
bound water layer and the diffuse bulk liquid, for 2 ; > sA. 
The barrier heights for the 100 and 111 surfaces are both 
large, at 9 and 10, respectively. These barriers are much 
larger then that depicted in Eig.j^a). This is a reflection 
of the reduced force needed for collective fluctuations to 
reorganize around the tagged molecule, rather than the 
forces to move a molecule relative to a static environ- 
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ment. 

We define indicator functions analogously as in the 
main text but now as a function of 2 ;, 

/li = 1 if < O.SA (Bl) 

= 0 else 

and 

gi = l a Zi> 3 A (B2) 

= 0 else 

and from these indicator functions we similarly compute 
a flux correlation function 

C(t) = , (b3) 

which is plotted in the inset to Fig.j^for both surfaces. In 
harmony with the result in the main text the asymptotic 
behavior of C{t) is kdt^ where kd is the rate constant 
using these basin definitions. As before this is exacted 
by a best fit line and results in kd = 0.053 ns“^ for the 
100 surface and kd = 0.028ns“^ for the 111 surface. 

As alluded to in the Section while use hi and gi 
rather then hi and gi does not change the overall rate 
constant, it does change the time for reaching steady- 
state. For the 100 surface this is significant, changing the 
lag time from it from 3ps to 20 ps. For the 111 surface 
it is not as significant, changing the mean time from 3ps 
to 12ps. This behavior is evident in the inset to Fig. 
where the linear growth regime can be gauged by the line 
of unit slope also plotted. The origin of this difference 
between the two surfaces is due to the 100 beiM signif¬ 
icantly more hydrophobic than the 111 surface,^ which 
thus enables much larger inter facial density fluctuations. 

Appendix C: Platinum 111 surface 

The main text focuses on the kinetics of the 100 sur¬ 
face, which is four fold-coordinated and more commensu¬ 
rate with extended hydrogen boning patterns resulting in 
unit monolayer coverage. We have also carried out this 
analysis for the 111 platinum surface, which is six-fold 
coordinated, and less commensurate with extended hy¬ 
drogen bonding patterns resulting in an average coverage, 
which is nearly 80%. Previous work has detailed differ¬ 
ences in the monolayer dynamics and solvation properties 
between these two surfaces. Here, we focus on only those 
differences that alter pathways for desorption. 

First, a number of quantitative differences result when 
desorption occurs at the 111 crystal facet compared to 
the 100 facet. As shown in Fig. the timescales are 
generically shorter with a rate 2 times larger and a mean 
instanton time 2.5 ps. These results are consistent with 
previous work that found faster relaxation times of wa¬ 
ter molecules adsorbed on the 111 compared to the 100. 
This is due to the six-fold coordination of the lattice with 




FIG. 9. (a) Probability distribution of surface reorganization 
energies. A, for the 111 surface. The dashed line is a Gaussian 
fit. (b) Probability distributions of F, in red, and G, in blue 
for the 111 surface. Dashed lines are Gaussian fits with same 
means and variances. 


is less commensurate with hydrogen bonding that intro¬ 
duces a larger degree of disorder on the 111 surface and 
facilities motion. 

These smaller characteristic timescales are mirrored in 
the relevant work quantities elucidated in the main text. 
Specifically, Fig. shows the distributions of A whose 
mean is similar to the 100 surface, but whose width is 
1.25 times larger leading to a high probability of lower 
values. Distributions of F and G are plotted in Fig. Ep- 
The mean values of each for the 111 surface are shifted 
relative to the 100 by about 2 /cbT. Within the limited 
statistics, these distributions are Gaussian. 

For the 111 surface, we And that the number of vacan¬ 
cies adjacent to a tagged site, increases monotonically 
with increasing reorganization energy. This is a reflec¬ 
tion of the stable hexagonal ring structure that water 
forms on the surface, which is commensurate with the 
6-fold coordination of the exposed lattice. For higher lo¬ 
cal surface densities, where tagged molecules do not sit 
adjacent to vacant sites, hydrogen bonding is frustrated 
and consequently these sites have a lower A. 
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